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Background 


Background 


• Modern simulation software is complex: 

• Implicit numerical methods 

• Massively parallel computers 

• Adaptive methods 

• Multiple, coupled physical processes 

• There are a host of existing software libraries that excel at treating 
various aspects of this complexity. 

• Leveraging existing software whenever possible is the most efficient 
way to manage this complexity. 
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Background 


Background 


• Modern simulation software is multidisciplinary: 

• Physical Sciences 

• Engineering 

• Computer Science 

• Applied Mathematics 


• It is not reasonable to expect a single person to have all the 
necessary skills for developing & implementing high-performance 
numerical algorithms on modern computing architectures. 

• Teaming is a prerequisite for success. 
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Background 


• A large class of problems are amenable to mesh based simulation 
techniques. 

• Consider some of the major components such a simulation: 
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Background 

• A large class of problems are amenable to mesh based simulation 
techniques. 

• Consider some of the major components such a simulation: 

O Read the mesh from file 
o Initialize data structures 

0 Construct a discrete representation of the governing equations 
O Solve the discrete system 
@ Write out results 

0 Optionally estimate error, refine the mesh, and repeat 
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• With the exception of step 3, the rest is independent of the class of 
problems being solved. 
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Background 

• A large class of problems are amenable to mesh based simulation 
techniques. 

• Consider some of the major components such a simulation: 

O Read the mesh from file 
o Initialize data structures 

0 Construct a discrete representation of the governing equations 
O Solve the discrete system 
@ Write out results 

0 Optionally estimate error, refine the mesh, and repeat 

• With the exception of step 3, the rest is independent of the class of 
problems being solved. 

• This allows the major components of such a simulation to be 
abstracted & implemented in a reusable software library. 
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The libMesh Software Library 


• In 2002, the libMesh library began with these ideas in mind. 

• Primary goal is to provide data structures and algorithms that can be 
shared by disparate physical applications, that may need some 
combination of 

• Implicit numerical methods 

• Adaptive mesh refinement techniques 

• Parallel computing 

• Unifying theme: mesh-based simulation of partial differential 
equations (PDEs). 
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The libMesh Software Library 


Key Point 

• The libMesh library is designed to be used by studendts, 
researchers, scientists, and engineers as a tool for developing 
simulation codes or as a tool for rapidily implementing a numerical 
method. 

• libMesh is not an application code. 

• It does not “solve problem XYZ.” 

• It can be used to help you develop an application to solve problem 
XYZ, and to do so quickly with advanced numerical algorithms on 
high-performance computing platforms. 
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Introduction 


Software Reusability 


• At the inception of libMesh in 2002, there were many high-quality 
software libraries that implemented some aspect of the end-to-end 
PDE simulation process: 

• Parallel linear algebra 

• Partitioning algorithms for domain decomposition 

• Visualization formats 

• ... 

• A design goal of libMesh has always been to provide flexibile & 
extensible interfaces to existing software whenever possible. 

• We implement the “glue” to these pieces, as well as what we viewed 
as the missing infrastructure: 

• Flexibile data structures for the discretization of spatial domains and 
systems of PDEs posed on these domains. 
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Introduction 


Software Reusability 



Library Structure 

• Basic libraries are 
libMesh’s “roots” 

• Application 
“branches” built off 
the library “trunk” 
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Library Design 


The “Glue 


• The C++ programming language provides a powerful abstraction 
mechanism for separating a software interface from its 
implementation. 

• The notion of Base Classes defining an abstract interface and Derived 
Classes implementing the interface is key to this programming model. 
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Introduction 


Library Design 


• The C++ programming language provides a powerful abstraction 
mechanism for separating a software interface from its 
implementation. 

• The notion of Base Classes defining an abstract interface and Derived 
Classes implementing the interface is key to this programming model. 

• The classic C++ example: Shapes. 

int main () 

{ 

/* using shapes polymorphically */ 

Shape * shapes [2]; 

shapes [0] = new Rectangle (10, 20, 5, 6); 
shapes [1] = new Circle (15, 25, 8); 

for (int i=0; i<2; ++i) 
shapes [ i ] ->Draw ( ) ; 
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Library Design 


Abstract Shape 


/* abstract interface declaration */ 
class Shape 
{ 

public : 

virtual void Draw () =0; 

virtual void MoveTo (int newx, int newy) = 0 
virtual void RMoveTo (int dx, int dy) = 0; 
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Library 


Specific Shape: Rectangle 


/* Class Rectangle */ 
class Rectangle : public Shape 
{ 

public : 

Rectangle (int x, int y, int w, int h) ; 
virtual void Draw (); 

virtual void MoveTo (int newx, int newy 
virtual void RMoveTo (int dx, int dy) ; 
virtual void SetWidth (int newWidth) ; 
virtual void SetHeight (int newHeight) ; 

private : 
int x, y; 
int width; 
int height; 
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Specific Shape: Circle 


/* Class Circle */ 
class Circle : public Shape 
{ 

public : 

Circle (int initx, int inity, int initr) 
virtual void Draw (); 

virtual void MoveTo (int newx, int newy) 
virtual void RMoveTo (int dx, int dy) ; 
virtual void SetRadius (int newRadius) ; 

private : 
int x, y; 
int radius; 
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Introduction 


Library Design 


int main () 

{ 

/* using shapes polymorphically */ 

Shape * shapes [2]; 

shapes [0] = new Rectangle (10, 20, 5, 6); 
shapes [1] = new Circle (15, 25, 8); 

for (int i=0; i<2; ++i) 

{ 

shapes [ i ] ->Draw ( ) ; 
DoSomethingWithShape (shapes [i] ) ; 

} 


/* access a rectangle specific function */ 
Rectangle rect(0, 0, 15, 15); 
rect . SetWidth (30); 
rect . Draw ( ) ; 



Library Design 


Examples of Polymorphism in 

4 


libMesh 
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Library Design 


The “Glue:” Linear Algebra 
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Library Design 


The “Glue:” I/O formats 
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Disretization: The Mesh 


libMesh::MeshBase 


libMesh::UnstructuredMesh 


z 


1 

1 


libMesh::ParallelMesh 


libMesh::SerialMesh 
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libMesh::Mesh 


libMesh: :BoundaryMesh 



Library Design 


Disretization: Geometric Elements 
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Disretization: Geometric Elements 
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Disretization: Geometric Elements 

















Library Design 


Disretization: Finite Elements 
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Library Design 


Algorithms: Domain Partitioning 
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Library Design 


Algorithms: Domain Partitioning 
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Library Design 


Algorithms: Error Estimation 
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Application Results 


Results from Physics Applications built 
on top of libMesh 
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Compressible Navier-Stokes 
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mass flux, H„ 


Motivation 


Application Results 


Arcjet Nozzle Calculation 
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Arcjet Nozzle Calculation 
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Coupled Pyrolysis, Temperature 
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Application Results 


Coupled Pyrolysis, Temperature 

Arc-.Tet 4" Iso-0 
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Coupled Pyrolysis, Pyrolysis gas mass flux, m 
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Coupled Pyrolysis, Pyrolysis gas mass flux, m 
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Application Results 


Coupled Thermal/Solid Mechanics 
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Coupled Thermal/Solid Mechanics 
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The MOOSE Framework - Gaston et al., INL 



MOOSE 

Derek Gaston, Idaho National Laboratory 

Collaborators: 

Dana Knoll, LANL 
Glen Hansen, INL 
Chris Newman, LANL 
Ryosuke Park, INL 
Cody Permann, INL 
David Andrs, INL 
Hai Huang, INL 
Michael Tonks, INL 
Robert Podgorney, INL 
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The MOOSE Framework - Gaston et al., INL 



MOOSE - Multiphysics Object Oriented 
Simulation Environment 

• A framework for solving computational nuclear engineering 
problems in a well planned, managed, and coordinated way 

- Leveraged across multiple programs 

• Designed to significantly reduce the expense and time 
required to develop new applications 

• Designed to develop analysis tools 

- Uses very robust solution methods 

- Designed to be easily extended and maintained 

- Efficient on both a few and many processors 

• Currently supports ~7 applications which are developed and 
used by ~20 scientists. 
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Application Results 


The MOOSE Framework - Gaston et al., INL 


MOOSE Ecosystem ^UU^Idaho National Laboratory 


Application 

Physics 

Start 

Time To 
Results 

Lines of Coc 


Thermo-mechanics, 




BISON 

Chemical Diffusion, 
coupled mesoscale 

June 2008 

4 Months 

931 

PRONGHORN 

Neutronics, Porous 
Flow, Eigenvalue 

September 

2008 

3 Months 

2,883 

SALMON 

Multiphase Porous 
Flow 

June 2009 

3 Months 

800 

MARMOT 

4 th Order Phasefield 
Mesoscale 

August 2009 

1 Month 

838 

RAT 

Porous ReActive 
Transport 

August 2009 

1 Month 

439 

FALCON 

Geo-mechanics, 
coupled mesoscale 

September 

2009 

3 Months 

810 


June 17, 2013 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite Element Library 



Application Results 


The MOOSE Framework - Gaston et al., INL 



BISON fuel performance 


• LWR, Triso, and TRU fuel performance code 

• Parallel 1D-3D thermomechanics code 

• Thermal, mechanical, and chemical models for FCI 

• Constituent redistribution 


• Material, fission product swelling, fission gas release models 

• Mesoscale-informed material models 
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Approach to Software Development 


Sidebar: 

Revision Control & Collaboration with GitHub 
Continuous Integration with Buildbot 
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A Generic Boundary Value Problem 


Solving Problems the libMesh way 
Discretizing a Generic Boundary Value Problem 
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A Generic BVP 


• We assume there is a Boundary 
Value Problem of the form 


du 

= F{u) 

G SI 

G{u) 

= 0 

G Q 

u 

= u D 

G 

N{u) 

= 0 

G 

!(JC,0) 

= uo(x) 




-w- 
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A Generic Boundary Value Problem 


• Associated to the problem 
domain Q is a libMesh data 
structure called a Mesh 

• A Mesh is essentially a 
collection of finite elements 

ri'-.= \Jn e 
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A Generic BVP 


• Associated to the problem 
domain Q is a libMesh data 
structure called a Mesh 

• A Mesh is essentially a 
collection of finite elements 

ri l := \Jn e 



• libMesh provides some simple structured mesh generation 
routines, file inputs, and interfaces to Triangle and TetGen. 
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Key Data Structures 


The Mesh Class 


The Mesh 


int main (int argc, char** argv) 

{ 

// Initialize the library. This is necessary because the library 
// may depend on a number of other libraries (i.e. MPI and PETSc) 
// that require initialization before use. When the LibMeshlnit 
// object goes out of scope, other libraries and resources are 
// finalized. 

LibMeshlnit init (argc, argv) ; 

// Create a mesh 
Mesh mesh; 

// Read the input mesh, 
mesh. read ("in.exo") ; 

// Print information about the mesh to the screen, 
mesh . print_inf o ( ) ; 

// Write the output mesh 
mesh. write ( "out . dat " ) ; 


} 
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Key Data Structures 


The Mesh Class 


The Mesh 


'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

* Running Example introduction_exl : 

* example-opt -d 3 reference_elements/3D/one_hex27 . xda -o output. xda 

•k-k-k-k-k-k-k-k-k-k'k'k'k'k'k'k'k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k'k'k'k'k'k-k-k-k-k-k-k-k'k'k'k-k-k'k-k-k-k-k-k-k-k 

Mesh Information: 
mesh_dimension ( ) =3 
spatial_dimension ( ) =3 
n_nodes ( ) =27 

n_local_nodes ( ) =27 
n_elem ( ) =1 

n_local_elem ( ) =1 
n_active_elem ( ) =1 
n_subdomains ( ) =1 
n_partitions ( ) =1 
n_processors ( ) =1 
n_threads ( ) =1 
processor_id ( ) =0 
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The Mesh Class 


Mesh Iterators 

void foo (const MeshBase &mesh) 

{ 

// Now we will loop over all the elements in the mesh that 
// live on the local processor. We will compute the element 
// matrix and right-hand-side contribution. Since the mesh 
// may be refined we want to only consider the ACTIVE elements, 
// hence we use a variant of the \p active_elem_iterator . 
MeshBase : : const_element_iterator 

el = mesh . active_local_elements_begin ( ) ; 

const MeshBase: : const_element_iterator 

end_el = mesh . active_local_elements_end ( ) ; 

for ( ; el != end_el; ++el) 

{ 

// Store a pointer to the element we are currently 
// working on. This allows for nicer syntax later, 
const Elem* elem = *el; 


June 17, 2013 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite Element Library 




The Mesh Class 


Mesh Iterators 


void foo (const MeshBase &mesh) 


// We will now loop over all nodes. 

MeshBase :: const_node_iterator node_it = mesh . nodes_begin ( ) ; 


const MeshBase :: const_node_iterator node_end = mesh . nodes_end ( ) ; 

for ( ; node_it != node_end; ++node_it) 

{ 

// the current node pointer 
const Node* node = *node_it; 
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The EquationSystems Class 


II Create a mesh. 

Mesh mesh; 

II Create a uniform 5x5 mesh on the unit square. 

MeshTools: : Generation : :build_square (mesh, 5, 5); 
mesh .print_inf o ( ) ; 

// Create an equation systems object. This object can contain 
II multiple systems of different flavors for solving loosely coupled 
II physics. Each system can contain multiple variables of different 
II approximation orders. The EquationSystems object needs a 
II reference to the mesh object, so the order of construction here 
// is important. 

EquationSystems equation_systems (mesh) ; 

// Now we declare the system and its variables. We begin by adding 
// a "TransientLinearlmplicitSystem" to the EquationSystems object, 

II and we give it the name "Simple System". 

equation_systems . add_system<TransientLinearImplicitSystem> ("Simple System") ; 

II Adds the variable "u" to "Simple System". "u" will be 
// approximated using first-order approximation. 

equation_systems . get_system ( " Simple System") . add_variable ("u", FIRST); 

// Next we'll by add an "ExplicitSystem” to the EquationSystems 
II object, and we give it the name "Complex System". 
equation_systems . add_system<ExplicitSystem> ("Complex System") ; 

II Give "Complex System" three variables — each with a different 
// approximation order. Variables "c" and "T" will use first-order 
// Lagrange approximation, while variable "dv" will use a 
II second-order discontinuous approximation space. 

equation_systems . get_system ( "Complex System") . add_variable ( "c" , FIRST); 
equation_systems . get_system ( "Complex System") .add_variable ("T", FIRST); 
equation_systems . get_system ( "Complex System" ). add_variable ( "dv" , SECOND, MONOMIAL); 

II Initialize the data structures for the equation system. 
equation_systems . init ( ) ; 

II Prints information about the system to the screen. 
equation_systems . print_inf o ( ) ; 
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a Background 

a The libMesh Software Library 
a Software Reusability 
a Library Design 
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a The Mesh Class 
a The EquationSystems Class 

$ Weighted Residuals 

^ Other Examples 



Weighted Residuals 


• The point of departure in any FE analysis which uses libMesh is the 
weighted residual statement 

(F(u),v)=0 VveV 


-w- 
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Weighted Residuals 


• The point of departure in any FE analysis which uses libMesh is the 
weighted residual statement 

(F(u),v)=0 VvgV 


• Or, more precisely, the weighted residual statement associated with 
the finite-dimensional space V h c V 


(F(u h ),v h ) = 0 Vv h £V h 
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Weighted Residuals 


Some Examples 


Poisson Equation 


—Au =/ g n 


June 17, 2013 50/73 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite Element Library 



Weighted Residuals 


Some Examples 


Poisson Equation 


— An =f g n 


Weighted Residual Statement 


(. F(m), v ) := / [Vw • Vv —fv\ dx 
J n 

+ / (' Vu ■ n) v ds 

J 
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Weighted Residuals 


Some Examples 


Linear Convection-Diffusion 

-kAu + b-X / u=f € L2 

- 
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Weighted Residuals 


Some Examples 


Linear Convection-Diffusion 

-kAu + b-X7u=f € L2 

Weighted Residual Statement 


(F(u),v) := / [kVu ■ Vv + (b ■ Vw)v — fv\ dx 

Jn 

+ k (Vm • n) v ds 
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Weighted Residuals 


Some Examples 


Stokes Flow 


Vp — vAu = / 

V • u = 0 


e Q 
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Weighted Residuals 


Some Examples 


Stokes Flow 


Vp — uAu 
V • u 


f 

0 


e Q 


Weighted Residual Statement 


u:=[u,p\ , v\=[v,q] 


(F(u),v) := / [— p (V • v) + uVu :Vv —f ■ v 

Jn 

+ {W -u)q\dx+ / (uVu — pi) n ■ v ds 

J d^N 
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Weighted Residuals 


Some Examples 


• To obtain the approximate problem, we simply replace u 4- u h , 
v 4- v h , and Q 4- Q h in the weighted residual statement. 
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Introduction 

a Background 

a The libMesh Software Library 
a Software Reusability 
a Library Design 

« Application Results 

a The Mesh Class 
a The EquationSystems Class 

^ Poisson Equation 



Poisson Equation 


Weighted Residual Statement 


• For simplicity we start with the weighted residual statement arising 
from the Poisson equation, with OCIn = 0, 


(F(u h ),v h ) := 



[Vu h • Vv /! -fv h ] dx = 0 


Vv' 1 e V h 


CFlRab 
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Poisson Equation 


Element Integrals 


• The integral over Q h . . . 


0 = f [Vu h ■ Vv /! -fv h \ dx Vv /! e V h 
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Poisson Equation 


Element Integrals 


• The integral over Q h 
finite elements: 


is written as a sum of integrals over the N e 


0 = / 

[Vu h ■ Vv /! -> /! " 

dx 

W 




II 

[Vu h • Vv /! -jv h ] 

dx 

W 


e=l 


e V h 
G V h 
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An element integral will have 
contributions only from the 
global basis functions 
corresponding to its nodes. 

We call these local basis 
functions (/>,-, 0 < i < N s . 


Finite Element Basis Functions 



Element Library 
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Poisson Equation 


Finite Element Basis Functions 


• An element integral will have 
contributions only from the 
global basis functions 
corresponding to its nodes. 

• We call these local basis 
functions 0 < i < N s . 


N s 


v In. = 


22 c ‘ ( * 

1=1 

N s . 

dx = 22 Ci / 4>i dx 


i=i 
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Poisson Equation 


Element Matrix and Load Vector 


The element integrals . . . 



Poisson Equation 


Element Matrix and Load Vector 


The element integrals . . . 



—Jv h ] dx 


are written in terms of the local “</><” basis functions 



Poisson Equation 


Element Matrix and Load Vector 


• The element integrals . . . 



—Jv h ] dx 


• are written in terms of the local basis functions 

y Uj / V</>; • Vd>i dx — / f4>idx , i = l,...,N„ 

Ja e Jn e 

• This can be expressed naturally in matrix notation as 

K e U e — F e 
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Global Linear System 


• The entries of the element stiffness matrix are the integrals 

K\j := [ V(/>j • V<t>i dx 

J 
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• The entries of the element stiffness matrix are the integrals 

K\j := [ V(/>j • V<t>i dx 

J 


• While for the element right-hand side we have 

: [ fOi dx 

J Vie 
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Global Linear System 


• The entries of the element stiffness matrix are the integrals 

K e i} := [ V(/>j • V<t>i dx 

J 


• While for the element right-hand side we have 

: [ fOi dx 

Jn e 


• The element stiffness matrices and right-hand sides can be 
“assembled” to obtain the global system of equations 


KU = F 
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• The integrals are performed on a “reference” element Cl e 


J 





June 17, 2013 58/73 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite Element Library 



Poisson Equation 


Reference Element Map 


• The integrals are performed on a “reference” element Cl e 


J 




Poisson Equation 


Reference Element Map 


• The integrals are performed on a “reference” element Cl e 


J 





• Chain rule: V = J := 




V 4>j ■ V 4>i dx 



%<t>j ■ \J\ d £ 
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Element Quadrature 


• The integrals on the “reference” element are approximated via 
numerical quadrature. 
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• The integrals on the “reference” element are approximated via 
numerical quadrature. 

• The quadrature rule has N q points “£ 9 " and weights "w q ". 
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Element Quadrature 


• The integrals on the “reference” element are approximated via 
numerical quadrature. 

• The quadrature rule has N q points “£ 9 " and weights "w q ". 


f\ = [mm 

J Ci e 

~ 5 ^/ (^(^9) )<^i (^9) I- 7 (C<?) I 

9=1 
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Element Quadrature 


• The integrals on the “reference” element are approximated via 
numerical quadrature. 

• The quadrature rule has N q points “£ 9 " and weights "w q ". 


K\j = [ Vrfj • \J\d£ 

JCl e 
N g 

9=1 
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libMesh Quadrature Point Data 


• libMesh provides the following variables at each quadrature point q 


Code 

Math 

Description 

JxW [q] 

|7(£ 9 )K 

Jacobian times weight 

phi [i] [q] 

MCq) 

value of i th shape fn. 

dphi [ i ] [q] 


value of i th shape fn. gradient 

xyz [q] 

*(6?) 

location of in physical space 
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• The libMesh representation of the matrix and rhs assembly is 
similar to the mathematical statements. 


for (q=0; q<Nq; ++q) 
for (i=0; i<Ns; ++i) { 

Fe(i) += JxW [q] *f (xyz [q] ) *phi [i ] [q] ; 

for ( j =0 ; j<Ns; ++j) 

Ke ( i , j ) += JxW [q] * (dphi [ j ] [q] *dphi [ i ] [q] ) ; 

} 
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Matrix Assembly Loops 


• The libMesh representation of the matrix and rhs assembly is 
similar to the mathematical statements. 


for (q=0; q<Nq; ++q) 
for (i=0; i<Ns; ++i) { 

Fe(i) += JxW [q] *f (xyz [q] ) *phi [i ] [q] ; 
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«=i 
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• The libMesh representation of the matrix and rhs assembly is 
similar to the mathematical statements. 
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for (i=0; i<Ns; ++i) { 
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m q )\w q 
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• The libMesh representation of the matrix and rhs assembly is 
similar to the mathematical statements. 
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//We now define the matrix assembly function for the Poisson system 
// by computing the element matrices and right-hand sides. We are 
// omitting BCs for brevity, 
void assemble_poisson (EquationSystemsS es, 

const std: : strings system_name) 


// Get a constant reference to the mesh object, 
const MeshBaseS mesh = es . get_mesh ( ) ; 


// The dimension that we are running 

const unsigned int dim = mesh.mesh_dimension ( ) ; 

// Get a reference to the LinearlmplicitSystem we are solving 

LinearlmplicitSystems system = es . get_system<Linear!mplicitSystem> ( "Poisson") ; 


//A reference to the DofMap object for this system. The DofMap 
// object handles the index translation from node and element 
// numbers to degree of freedom numbers. We will talk more about 
/ / the DofMap in future examples . 
const DofMapS dof_map = system. get_dof_map(); 


// Get a constant reference to the Finite Element type for the first 
// (and only) variable in the system. 

FEType fe_type = dof_map . variable_type (0) ; 


// Build a Finite Element object of the specified type. Since the 
// FEBase : :build ( ) member dynamically creates memory we will store 
// the object as an AutoPtr<FEBase> . This can be thought of as a 
// pointer that will clean up after itself. 

AutoPtr<FEBase> fe (FEBase : :build (dim, fe_type) ) ; 

//A 5th order Gauss quadrature rule for numerical integration. 
QGauss qrule (dim, FIFTH) ; 
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Matrix Assembly Loops 


// Tell the finite element object to use our quadrature rule. 
fe->attach_quadrature_rule (Sqrule) ; 

// Here we define some references to cell-specific data that will be 
// used to assemble the linear system. We begin with the element 
// Jacobian * quadrature weight at each integration point, 
const std: :vector<Real>& JxW = fe->get_JxW ( ) ; 

// The physical XY locations of the quadrature points on the 
// element. These might be useful for evaluating spatially varying 
// material properties at the quadrature points, 
const std: :vector<Point>& q point = f e->get_xyz ( ) ; 

// The element shape functions evaluated at the quadrature points, 
const std: :vector<std: : vector<Real> >& phi = fe->get_phi ( ) ; 

// The element shape function gradients evaluated at the quadrature 
// points. 

const std: :vector<std: :vector<RealGradient> >& dphi = f e->get_dphi ( ) ; 

// Define data structures to contain the element matrix and 
// right-hand-side vector contribution. Following basic finite 
// element terminology we will denote these "Ke" and "Fe". 
DenseMatrix<Number> Ke; 

DenseVector<Number> Fe; 

// This vector will hold the degree of freedom indices for the 
// element. These define where in the global system the element 
// degrees of freedom get mapped, 
std: : vector<dof_id_type> dof_indices; 

// Now we will loop over all the elements in the mesh. We will 
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Matrix Assembly Loops 

// compute the element matrix and right-hand-side contribution. See 
// example 3 for a discussion of the element iterators. 

MeshBase : : const_element_iterator el = mesh . active_local_elements_begin ( ) ; 

const MeshBase :: const_element_iterator end_el = mesh . active_local_elements_end ( ) ; 

for ( ; el != end_el; ++el) 

{ 

// Store a pointer to the element we are currently working on. 

// This allows for nicer syntax later, 
const Elem* elem = *el; 

// Get the degree of freedom indices for the current element. 

// These define where in the global matrix and right-hand-side 
// this element will contribute to. 
dof_map . dof_indices (elem, dof_indices) ; 

// Compute the element-specific data for the current element. 

// This involves computing the location of the quadrature points 
// (q point) and the shape functions (phi, dphi) for the current 
// element. 
fe->reinit (elem) ; 

// Zero the element matrix and right-hand side before summing 
// them. We use the resize member here because the number of 
// degrees of freedom might have changed from the last element. 

// Note that this will be the case if the element type is 
// different (i.e. the last element was a triangle, now we are 
// on a quadrilateral) . 

Ke. resize (dof_indices . size ( ) , 
dof_indices . size ( ) ) ; 

Fe. resize (dof_indices . size ( ) ) ; 
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// Now we will build the element matrix. This involves a double 
// loop to integrate the test funcions (i) against the trial 
// functions (j) . 

for (unsigned int qp=0; qp<qrule . n_points ( ) ; qp++) 
for (unsigned int i=0; i<phi . size ( ) ; i++) 
for (unsigned int j=0; j<phi . size ( ) ; j++) 

Ke ( i , j ) += JxW[qp] * (dphi [i] [qp] *dphi [ j] [qp] ) ; 

// Now we build the element right-hand-side contribution. This 
// involves a single loop in which we integrate the "forcing 
// function" in the PDE against the test functions, 
for (unsigned int qp=0; qp<qrule . n_points ( ) ; qp++) 
for (unsigned int i=0; i<phi . size ( ) ; i++) 

Fe(i) += JxW [qp] *10 . *phi [i] [qp] ; 

// If we are using an adaptive mesh this will apply any hanging 
// node constraint equations 

dof_map . heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices) ; 

// The element matrix and right-hand-side are now built for this 
// element. Add them to the global matrix and right-hand-side 
// vector. The SparseMatrix : : add_matrix ( ) and 
// NumericVector : : add_vector ( ) members do this for us. 
system. matrix->add_matrix (Ke, dof_indices) ; 
system. rhs->add_vector (Fe, dof_indices) ; 

} 

} 


June 17, 2013 65/73 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite Element Library 



Other Examples 


a Background 

a The libMesh Software Library 
a Software Reusability 
a Library Design 

« Application Results 

a The Mesh Class 
a The EquationSystems Class 

Q Other Examples 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite 



CFIjLab 


Library 


June 17, 2013 


66/73 



Other Examples 


Convection-Diffusion Equation 




The matrix assembly routine for the linear convection-diffusion 
equation, 

—kAu + b ■ Vu =/ 


for (q=0; q<Nq; ++q) 

for (i=0; i<Ns; ++i) { 

Fe(i) += JxW[q] *f (xyz [q] ) *phi [i] [q] ; 


for (j = 0; 
Ke (i, j ) 


} 


j<Ns; ++j) 

+= JxW [q] * (k* (dphi [ j ] [q] *dphi [ i ] [q] ) 
+ (b*dphi [ j] [q] ) *phi [i] [q] ) ; 
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Convection-Diffusion Equation 




The matrix assembly routine for the linear convection-diffusion 
equation, 

—kAu + b ■ Vu =/ 


for (q=0; q<Nq; ++q) 

for (i=0; i<Ns; ++i) { 

Fe(i) += JxW[q] *f (xyz [q] ) *phi [i] [q] ; 


for (j = 0; 
Ke (i, j ) 


} 


j<Ns; ++j) 

+= JxW [q] * (k* (dphi [ j ] [q] *dphi [ i ] [q] ) 
+ (b*dphi [ j] [q] ) *phi [i] [q] ) ; 
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Convection-Diffusion Equation 




The matrix assembly routine for the linear convection-diffusion 
equation, 

—kAu + b • Vw =/ 


for (q=0; q<Nq; ++q) 

for (i=0; i<Ns; ++i) { 

Fe(i) += JxW[q] *f (xyz [q] ) *phi [i] [q] ; 


for (j = 0; 
Ke (i, j ) 


} 


j<Ns; ++j) 

+= JxW [q] * (k* (dphi [ j ] [q] *dphi [ i ] [q] ) 
+ (b*dphi [ j] [q] ) *phi [i] [q] ) ; 
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Stokes Flow 


• For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

U\ 

K e 

IS ~U2U\ 

K e 

n 'UiU2 

K e 

iY M2^2 

K e 

uip 

K lp 


r u e i 

U U 1 

u e 


r F e 

«i 

F e 

«2 

K p Ui 

Ku2 

K PP \ 


u e 

L p J 


F e 
L p 


• We have an array of submatrices: Ke [ ] [ ] 


4 
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Stokes Flow 


• For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

^ li\U\ 

K e 

I ^U2U\ 

K e 

I ^UiU2 

K e 

1 ^U2U2 

K e 

uip 

K lp 


r u e i 

U U 1 

u e 


r F e 

«i 

F e 

«2 

K P U1 

K PU 2 

Kp \ 


u e 

L p J 


F e 
L p 


• We have an array of submatrices: Ke [ 0 ] [ 0 ] 


4 
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For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

K e 

rS 'UiU2 

K e 

in p 

K e 

iV M2Wl 

K e 

Ke u,p 

K Un 


Kp 


U e 

U U 1 

u e 

U2 


ut 


We have an array of submatrices: Ke [ 1 ] [ 1 ] 


F e 

«i 

F e 

»2 
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Other Examples 


For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

K e 

n 'UiU2 

K e 

u\p 

K e 

I ^U2U\ 

K e 

lV U2U2 

K€ u,p 

K Un 

Ku 2 

Kr 


U e 

U U 1 

u e 

U7. 


ut 


We have an array of submatrices: Ke [2 ] [ 1 ] 


F e 

«i 

»2 


Kirk, Peterson, Stogner (NASA, INL, UT) The libMesh Finite Element Library 


CFlUab 

June 17,2013 68/73 
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For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

K e 

n 'UiU 2 

K e 

U\P 

K e 

I ^U 2 U\ 

K e 

iY M2^2 

K U,P 

Ktn 

Ktt 2 

Kp 


U e 

U U 1 

u e 

U7. 


Ui 


F e 

« i 

F e 

«2 


And an array of right-hand sides: Fe 
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For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

*^U\ U\ 

K e 

n 'UiU2 

K e 

Ml P 

K e 

I ^U2U\ 

K e 

iY M2^2 

K, P 

Km 

Km 

Kp 


U e 

U U 1 

u e 

U7. 


ut 


And an array of right-hand sides: Fe [ 0 ; 


F e 
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F e 
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Other Examples 


For multi-variable systems like Stokes flow, 


G 0 C 


Vp — z'Ah = / 

V • u = 0 

The element stiffness matrix concept can extended to include 
sub-matrices 


K e 

U\ 

K e 

n 'UiU2 

K e 

U\P 

K e 

I ^U2U\ 

K e 

iY «2^2 

K U,P 

Ktn 

Ktn 

Kr 


U e 

U U 1 

u e 

U7. 


ut 


And an array of right-hand sides: Fe [ 1 ; 


F e 

u i 
F e 
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Other Examples 


Stokes Flow 


• The matrix assembly can proceed in essentially the same way. 

• For the momentum equations: 


for (q=0; q<Nq; ++q) 
for (d=0; d<2; ++d) 

for (i=0; i<Ns; ++i) { 

Fe[d](i) += JxW[q] *f (xyz [q] , d) *phi [i] [q] ; 


} 


for (j=0; j<Ns; ++j) 

Ke [d] [d] (i, j) += 

JxW [q] *nu* (dphi [ j ] [q] *dphi [ i ] [q] ) ; 


J 


June 17,2013 69/73 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite Element Library 



Some Extensions 


Introduction 

a Background 

a The libMesh Software Library 
a Software Reusability 
a Library Design 

Q Motivation 

« Application Results 

a The Mesh Class 
a The EquationSystems Class 

^ Some Extensions 


Kirk, Peterson, Stogner (NASA, INL, UT) 


The libMesh Finite 


Library 


June 17, 2013 


70/73 



Some Extensions 


Time-Dependent Problems 


• For linear problems, we have already seen how the weighted residual 
statement leads directly to a sparse linear system of equations 


KU = F 
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Time-Dependent Problems 


For time-dependent problems, 


du CY \ 

a, = F{ “> 


we also need a way to advance the solution in time, e.g. a 0-method 


u n+l - u n 
At : 



ue 


( F(u 0 ),v h ) Vv /! G V h 
du' ,+l + (1 - 9) u n 


• Leads to KU = F at each timestep. 


CFfflLab 
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Nonlinear Problems 


• For nonlinear problems, typically a sequence of linear problems must 
be solved, e.g. for Newton’s method 

(F\u k )Su k+1 ,v) = -( F(u k ),v ) 

where F'{u k ) is the linearized (Jacobian) operator associated with the 
PDE. 

• Must solve KU = F (Inexact Newton method) at each iteration step. 
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